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Active particles under confinement 
1. Introduction 



Aggregation of motile organisms, or self-propelled particles in general, is a familiar 
phenomenon, which is currently driving an intense interest in many scientific disciplines 
that include physics, biology and sociology pp. Many different physical mechanisms 
can lead to the aggregation of active systems. For example, these mechanisms can 
be directional alignment interactions [2j |3j HJ El E], density-dependent motility [7J [8], 
and aggregation at the wall when the system is under confinement P, QUI ELD, fT2] . 
In this theoretical paper, I will focus on the latter mechanism: confinement-induced 
aggregation of active particles. Experimentally, this mechanism is of particular interest 
given its inevitability due to the finite size of an experimental apparatus. The advent of 
microfluidic devices as an experimental tool to study active systems further motivates 
the current study since the channel walls constitute an integral part of the system [13J. 
In the classical theory of liquids, aggregation of particles close to the wall is well 
understood and results from the pairwise interactions, such as volume exclusion, of 
the particles in the system (e.g., see [H] for a review). Density functional theory 
is an effective way of investigating such a phenomenon, and a parallel effort has 
been pursued recently for a system of self-propelled rods [10J. On the other hand, 
if pairwise interactions between the particles become negligible, such a confinement- 
induced aggregation phenomenon will cease to exist in an equilibrium system, while 
it will continue to persist in an active system. This is the effect I will focus on here. 
Specifically, I will investigate a system of active particles that travel at a fixed speed 
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Figure 1. (a) A schematic depiction of the model system considered in this paper. 
Active particles are confined in a channel that is infinite in the y-direction but bounded 
in the x-direction. The walls are assumed to be perfectly rigid in that the particles 
cannot penetrate them but the particles are free to rotate at and slide along the walls. 
The channel inside is allowed to have two distinct rotational diffusion constants: a\ 
for x > and 02 for x < 0. (b) The simplified model also considered here restrict the 
active particles to have only six distinct directions of travel. 



Active particles under confinement 3 

in a channel, undergo rotational diffusion, and react only with the walls. I will derive 
from first principles a set of equations that govern the density profile of this system, 
and show that the density function can be written as a series expansion of the Mathieu 
functions |15| . To make further analytical process, I will introduce a simplified model 
in which the particles' directions of travel are discretised. 

Besides determining quantitatively how confinement induces aggregation, I will also 
discuss the gradient formation of particle density in a channel that is partitioned into 
two regions within which the active particles exhibit distinct rotational fluctuations 
(Fig. Ufa)). Such a scenario can for instance be engineered by having two background 
mediums coexisting in a microfluidic channel, and has been recently investigated in a 
combined experimental and numerical work |16| . 

The structure of this paper is as follows: I will first specify the model under study 
in Section El and present three "predictions" based on intuitive arguments in Section El 
The mathematical formulation of the system will be discussed in Section HI A simplified 
model and its full analysis will be presented in Section [51 The paper ends with a 
discussion of the findings. 

2. Model 

I consider a minimal model of active particles in two dimensions, where every particle 
is assumed to have constant speed u. Noise, of strength a, is incorporated in the 
direction of travel and the particles interact only with the boundary of the system. 
Mathematically, the equations of motion for such a iV-particle system are 

^j = t*v(ft) (1) 

^ = yfcntf) (2) 

where 1 < % < N, R = (r 1 ,...,r N ), 6 = {8 1 ,...,6 N ), v(0) = (cos 6, sin 0), and (Eq. W) 
is an Ito stochastic differential equation [17] such that 

(*(*)> = , (Vi(t)rj j (t > )) = S ij S(t-t > ). (3) 

The specific geometry considered is depicted in Fig. [Ha) where the walls are assumed 
to be perfectly rigid and smooth such that the particles are free to rotate at and slide 
along the walls. Note that in this model, I ignored the thermal translational diffusion D 
of the particles [T8l [19] . As explained in [20], this simplification is valid if D <C v?a~ x , 
which applies, e.g., for the bacterium Escherichia coli. 

Without loss of generality, I will from now on set the length scale and temporal 
scale so that the speed u is one and the width of the channel 2d is 2, although I will 
continue to keep the symbols in the equations so that the dimensions of the terms in 
the equations are always clear. 
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3. Intuitive arguments 

I will first analyse the system based on physical arguments in this section. Since there are 
no interactions between particles by assumption, I can focus on the average behaviour 
of a single particle. Specifically, denoting the average time the average particle spent at 
the wall by t w , and the average time spend in the channel by r c , the boundary condition 
suggests that if a particle is stuck at a wall, i.e., having its direction of travel pointing 
towards the wall, then it will only leave the wall if it manages to rotate its direction to 
point away from it, so the average time spent at the wall can be estimated as the time 
it takes to rotate, say, by one radian. As a result, 

T w ~ l/o . (4) 

To estimate r c , one has to worry about whether the particles in the channel is dominantly 
ballistic or diffusive. These two behaviours are dictated by the rotational diffusion 
constant a. In other words, if a is large, then the particles will have rotated around 
many times over before reaching the walls and so the motion of the particle is diffusive 
with an effective diffusion constant -D e // = u 2 /a (see, e.g., Ch. 6 in [21J). On the other 
hand, if a is small, the motion of the particle inside the channel will be predominantly 
ballistic. Therefore, 

Tc „ J d2 l D zff ' diffusive 
| d/u , ballistic 

Remembering that the units are set such that d = u = 1, r c becomes 

r c ~j a ' a>>1 (6) 

Together with (Eq. HJ), one may estimate the ratio of the number of particles inside the 
channel n c over the number of particles at the walls n w to be 

a< 1 . U 

If the channel is now partitioned into two regions system with distinct rotational 
diffusion constants a\ and a 2 (Fig. Ufa)), an interesting question that is of experimental 
interests [16] is which side of the region will be of higher density. Focusing on the 
Oi, 02 3> 1 regime, one may expect that the particles would behave diffusively inside the 
channel, with effective diffusion constant D\ ~ 1/ai and D2 ~ 1/02 in the two distinct 
regions. Intuitively, if there is a gradient in the system, then the region with the lower 
diffusion constant would have a higher number of particles since the slow mobility would 
serve to localise particles (TJ [8] . Interestingly, such an effect has also been observed in the 
embryo of the round worm Caenorhabditis elegans, where regions with lower diffusion 
constants serve to localise certain diffusing proteins [22j [23]. If this expectation extends 
to our active system, then one prediction would be that for a%, ai ^> 1, a\ < a? implies 
n^p < n^ where n® denotes the number of particles in region i. 
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Summarising this section, the following three predictions seem to follow from the 
intuitive arguments presented above: 

(i) In the symmetric case (a\ — a 2 = a), if a -C 1, then n c /n w ~ a (?) 
(ii) In the symmetric case, if a 3> 1, then n c /n w ~ a 2 (?) 
(hi) In the asymmetric case, if ai < a 2 and ai, 02 3> 1, then n^ < n^ (?) 

With the mathematical tools developed in the next section, I will show that in fact only 
the first "prediction" is true. 

4. Mathematical formulation 

I will now provide a mathematically rigorous formulation of the problem. If we denote 
the probability distribution of the density of particles in the state (R, 0) at time t by 
f(t, R, O), then the Fokker-Planck equation corresponding to the system is [24J: 



f = E{«J/-«V r ,.[v( W ]}. (8) 

Since by assumption, the particles do not interact with each other, the dynamics of the 
system is fully captured by the corresponding single-particle density function p(ri, 9 2 ) = 
N J dr 2 ■ • • r N d9 2 ■ ■ ■ 9 N f(Il, 6), which is governed by the following dynamical equation 
[2S1 E5J: 

% = a W- mV - • ( v (^ • ( 9 ) 

Focusing on the steady-state solution and assuming the channel is infinite along the 
y-direction, the equation simplifies to 

d 2 p a dp 

a W = uco * e T x - (10) 

Employing the separation of variables method by assuming that p(x,6) = A(6)B(x), I 
arrive at the two equations 

d ee A = q cos 9 A , (11) 

d x B = qaB , (12) 

where q corresponds to a set of constants to be determined. Eq. f|T2|) can be easily 
solved: 

B = Ke Dqx , (13) 

while Eq. ( ITTj) indicates that the function A is a Mathieu function, usually denoted as 
ce2n(#, <?2n) where n G N, such that the corresponding characteristic number is zero [13] . 
Since the characteristic numbers for these types of Mathieu functions are identical for 
±g, the expansion for p is 

p(x, 6) = J2 [a n e q2 " ax ce 2n (9/2, 2q 2n ) + b n e-^ ax ce 2n (9/2, -2q 2n )] . (14) 

n>0 
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Figure 2. Steady-state distribution of active particles in the symmetric case (a\ — 02). 
(a) The density map of the active particles in the system at the steady state as a 
function of position and direction based on extensive Brownian dynamics simulation 
[2 7) . The highest density is arbitrarily set to one (see colorbar). (b) The normalised 
distribution of the particle's direction at the left wall L(9). (c) The normalised density 
of particles within the channel Ptot{%) = J d9p(x, 9). (d) The ratio of the number of 
particles in the channel vs. that at the wall (n c /n w ) as a function of a, for both the 
full and simplified 6-direction models, the curve depicts the analytical expression in 
Eq. (HU). 



As usual, the expansion parameters are obtained by the boundary conditions, which I 
will determine now. 

Denoting the number of particles at the wall on the left (x = —d) travelling in 
the direction 9 by L(0), L(9) must be zero for — tt/2 < 9 < tt/2 as particles travelling 
in those directions would not stay at the wall (Fig. [21(b)). For rr/2 < 9 < 3rr/2, the 
equation for L(9) can be derived by the conservation of particles as follows: The flux 
of particles that arrive at the wall is —ucos9p L (9) where p L (9) denotes p(x = —d,9). 
At the wall, the particles that are stuck (i.e., particles with direction —tt/2 < 9 < tt/2) 
undergo rotational diffusion. Therefore, the density function L(9) satisfies the following 
differential equation: 

d 2 L 

a —- = ucos9p L {9) , tt/2 < 9 < 3tt/2 , (15) 

o9 l 
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with the boundary conditions L{9 = ±7i"/2) = 0. Finally, the flux of particles that leave 
the wall corresponds to the flux of particles that reach the angle 9 = ±7r/2. For instance, 
for the out-flux source at 9 = tt/2, the magnitude of the out-flux is \aOL{9 = tt/2)/09\. 
Together with the out-flux source at 9 = — tt/2, the overall out-flux at the left wall is 



a<6(6-ir/2) 



0L(9 = tt/2) 



09 



+ 6(9 + tt/2) 



dL(9 = -tt/2) 



09 



(16) 



Similarly, the equation governing the steady-state profile at the wall on the right 
R(9) is 

"j R (9) , -tt/2 < 9 < tt/2 , (17) 



09 2 



-U COS I 



with the same boundary conditions: R(9 = +tt/2) = 

As aforementioned, at the steady-state, the particles that are accumulated at the 
walls eventually re-enter the channel at four delta sources (two at either wall) (see Fig. 
[2]^a)). To account for these influxes of particles from the walls, Eq. (fit)]) is modified as 
follows: 



Op Op 

i—— = ucos9— — \-aS(x— d)5(9±Tr/2) 
09 2 Ox 



OR 



09 



+aS(x+d)5(9±Tr/2) 



OL 



09 



.(18) 



In summary, the steady-state density profile of the active particles is obtained by solving 
the coupled differential equations Eqs (!T5|) . ( 1T7|) and ( TT8"]) . 

The standard way to proceed at this point is to substitute into the above coupled 
differential equations the expansion for p in terms of the Mathieu functions (Eq. 
( Tl4l) ). and then try to ascertain the expansion coefficients accordingly. In practice, 
further analytical process is difficult since manipulations of the Mathieu functions are 
analytically intractable. Indeed, the numerical challenge to determine the characteristic 
numbers for Mathieu functions is comparable to performing Brownian dynamics 
simulation of the system [T5], |27] . Therefore, in order to further reveal the underlying 
physics of this problem, I will now introduce and analyse a simplified version of the 
original model that admits full analytical solution. 



5. Simplified 6-direction model 

In this simplified model, the angular components are partitioned into six directions as 
shown in Fig. [D(b) and each particle will switch to a neighbouring direction of travel 
with rate a. I will denote the density function inside the channel travelling towards the 
z-th direction by Pi, and those at the wall on the right and on the left by Ri and L t 
respectively. The translational symmetry along the y direction indicates that p± = P3, 
Pe — Vii an d also that the functions p only depend on x but not y. As a result, the 
dynamical equations for the system are 

OtPi = a(p 2 +P6- 2pi) - bO x pi (19) 

t p 2 = 2a(pi - p 2 ) - x p 2 (20) 

t p 5 = 2a(p & - p 5 ) + x p 5 (21) 
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d t p 6 = a(p 1 +p 5 - 2p 6 ) + bd x p 6 
a/3/ 2 since b = cos it/3. In matrix form, dp/dx 



8 

(22) 



where b = 
(pi,P2,P5,P6) T and 



Mp where p 



M 



( -2a /b a/b 

2a -2a 



V -a/b 



(23) 



a/b \ 



2a -2a 

-a/b 2a/b J 

The four eigenvalues associated to the matrix M are {0, 0, —07/6, 07/6} where 7 = 
V3 + 46 + 46 2 , and the corresponding eigenvectors are the column vectors in the 
following matrix: 

2 + 7 -^ 1 \ 

-4(6 + b 2 ) 



V 



/o 1 

1 
1 

V° 1 



26 

26+7 

1 



26+7 
_ 26(2-7) 

26+7 



26 

26+7 



_ 4(& + 6 2 ) _ 26(2-7) 



26+7 



2 + 7 



26 

26+7 



(24) 



/ 



Furthermore, the steady-state solution is of the form: 

Ph{x) 



Y. c ^ kX v hk 



(25) 



where o~k G {0, 0, —07/6, 07/6} are the corresponding eigenvalues of M. The above 
expression is to be compared to Eq. (fT4"j) in the continuum model. 

To determine the coefficients Ck in Eq. ( 1251) . I employ the discrete version of the 
boundary conditions as discussed in the previous section (Eq. ( JT5"|) ). which for the wall 
on the left gives, 






L 2 
u 

a 







and the flux of particles out of the wall on the left is (see Eq. ([16 
aL 6 



v\ 



bu 

Yb^ 



(26) 

(27) 

(28) 



(29) 
(30) 



where the second equality follows form Eq. (|28|) . Similar expressions can be readily 
found for the Ri. In particular, the condition in Eq. (|30|) and the corresponding one for 
pf enable me to determine the coefficients in the expansion in Eq. ( )25l) completely. 



5.1. Symmetric case (a\ = a 2 = a) 

In the symmetric case, the interesting quantity to investigate is the ratio of the number of 
particles inside the channel n c vs. the number of particles at the wall n w . The analytical 



Active particles under confinement 

(a) 








(c) 




0.56 










i 


+ Simulation 




0.54 


v 


Theory 




2 0.52 


\ 






o 


\ I 


0.5 


V / 


0.48 


w 



Figure 3. Steady-state properties of the simplified system in the asymmetric case 
where two regions with distinct rotational diffusion constants co-exist. The results are 
based on the analytical theory presented in the main text, (a) The contour plot of 
the ratio of number of particles inside the left-half of the channel vs. the right-half 
of the channel excluding the particles at both walls ( n c /rv c M as a function of a\ 
and a-i. (b) The contour plot of the ratio of number of particles inside the left-half 
of the channel vs. the right-half of the channel including the particles at both walls 
( (L tot + n c )/(Rtot + n c ) ) as a function of a\ and a-i. (c) The gradient of particle 
density (ptot(x) — ^2kPk(%)) formed inside the channel for the case a\ = 0.2 and 
a 2 = 5. 



expression for n c /n w can be obtained straight-forwardly with the help of Mathematica 
[28] . Although the analytical expression is complicated and difficult to decipher, it does 
provide the following asymptotic expressions with respect to a: 

1.54a , a < 1 



n„ 



3b a ~ 

12a(-l+7+2b7+7 2 



1.17a 



a> 1 



(31) 



7(7+126 2 +4 7 +2b(7+37)) 

This result indicates that the ratio increases linearly in both asymptotic regimes, albeit 
with two different slopes. These analytical results are in perfect agreement with results 
from Brownian dynamics simulation (Fig. El^d)). Referring to Section El expectation (i) 
is thus verified while expectation (ii) is invalidated. I will come back to discuss why this 
is the case in the Discussion Section. 



5.2. Asymmetric case (a\ ^ 02) 

In the two- region case as depicted in Fig. Ufa), the expansion coefficients in Eq. (1251) 
are obtained by employing the boundary conditions at both walls (Eq. (1301)). and by 
matching the densities Pi at x = 0. Note that the matching condition at x = ensures 
the continuity of the density functions but not their slopes. This is fully supported by 
simulation result (Fig. Etc)). 

Fig. Etc) also shows that there is a gradient of particle density inside the channel, 
one natural question is which region possesses more particles. The answer can again be 
obtained from the analytical solution for pi, i?, and Lj, and the result is that the region 
with a larger rotational diffusion constant always has a smaller number of particles, 
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whether the comparison concerns the particles inside the channel only (Fig. |3Ja)) or 
includes the particles at the walls (Fig. [3(b)). In other words, if a\ < a 2 , then n^' > n^ . 
Hence, the expectation (iii) discussed in Sect. |3jis not valid. Note that this surprising 
phenomenon has also been observed previously in a simulation study of a similar type 
of active system [16] . 

6. Discussion 

As we have seen, among the three expectations arisen from intuitive arguments in Sect.[3j 
only the first expectation is verified. The reason that expectation (ii) is invalid is due 
to the overall conservation in the number of particles in the system. In other words, 
although it remains true that as a — > oo, n w ~ 1/a, n c cannot scale with a as this would 
lead to a divergence in the total number of particles. So as a increase, n c also increases 
but eventually saturates to a number bounded below by the total number of particles 
in the system. This is why n c /n w remains linear in a even in the large a regime. 

In the asymmetric case, in contradiction to the expectation (iii) from Sect. [3j a 
higher a<± in relation to a\ always leads to a higher density of particles in region 1, 
irrespective of whether ax, a^ ^> 1 or not (see Fig. [3ja)). A physical explanation may be 
that first, a higher a leads to a lower accumulation of particles at the wall (Eq. H]); and 
second, the exponent (i.e., the eigenvalues) scales with a (Eq. (1251) and so a higher a will 
lead to more rapid decay of density from the wall (Fig. [3jc)). These two complimentary 
effects turn out to be dominant in determining the gradient of the particle density, and 
thus lead to the observation that a higher a always leads to a lower density of particle 
in that region. Of course, the validity of this argument is only confirmed by previous 
analytical investigation. 

In summary, I have studied the confinement- induced aggregation phenomenon 
in a system of self-propelled particles inside a channel. The model system under 
consideration is drastically simplified in order to focus on the effects of confinement 
on aggregation. In particular, the particles are assumed to be non-interacting and the 
walls are assumed to be perfectly rigid and smooth. Even in this simplified scenario, 
solving the rotationally-continuous model is mathematically difficult and thus forbids 
an analytical understanding of the system. Therefore, I introduced a further-reduced 
model where the particles can only take six different directions of travel. This enabled 
me to solve the model analytically and to prove that when the rotational fluctuation 
a is uniform in the system, the ratio of the number of particles inside the channel vs. 
that accumulated at the walls increases linearly with a in both the aC 1 and a ^> 1 
regimes. I also demonstrated that if one side of the channel has a larger rotational 
fluctuation, this inevitably leads a lower density of particles on that side. Given the 
fundamental nature of the model and the general nature of the observations, the effects 
of confinement on active systems discussed here should continue to be present in more 
elaborate models and in real experimental systems. In particular, it would be interesting 
to study to what extend confinement- induced gradient formation in the bulk of an active 
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system can lead to the emergence of pattern formation in the whole system. 
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